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Rapidly rotating, slightly non-axisymmetric neutron stars emit nearly periodic gravitational waves 
(GWs), quite possibly at levels detectable by ground-based GW interferometers. We refer to these 
sources as "GW pulsars" . For any given sky position and frequency evolution, the J^-statistic is the 
optimal (frequentist) statistic for the detection of GW pulsars. However, in an "all-sky" searches 
for previously unknown GW pulsars, it would be computationally intractable to calculate the (fully 
coherent) ^-"-statistic at every point of (suitably fine) grid covering the parameter space: the number 
of gridpoints is many orders of magnitude too large for that. Therefore, in practice some non-optimal 
detection statistic is used for all-sky searches. Here we introduce a "phase-relaxed" J-"pr-statistic, 
which we denote J-pr, for incoherently combining the results of fully coherent searches over short 
time intervals. We estimate (very roughly) that for realistic searches, our J^pr is ~ 10 — 15% more 
sensitive than the "semi-coherent" J^-statistic that is currently used. Moreover, as a byproduct of 
computing J-pr, one obtains a rough determination of the time-evolving phase offset between one's 
template and the true signal imbedded in the detector noise. Almost all the ingredients that go 
into calculating J^pr are already implemented in the LIGO Algorithm Library, so we expect that 
relatively little additional effort would be required to develop a search code that uses J-pr- 

PACS numbers: 95.55.Ym, 04.80.Nn, 95.75.Pq, 97.60. Gb 



I. INTRODUCTION 

The J^-statistic is the optimal frequentist statistic 
for the detection of nearly monochromatic gravitational 
waves (GWs) from a neutron star with known (or as- 
sumed) sky location and frequency evolution [ij. The 
basic idea behind the J^-statistic is simply this: for any 
given sky location and frequency evolution, the set of pos- 
sible GW signals forms a four-dimcnsional (real) vector 
space P, Q ■ The four basis vectors arc the two quadra- 
tures (sin and cos) of each of the two polarization bases, 
4- and x . Because the set is a vector space (not just a 
4-d manifold), it is computationally trivial to maximize 
the likelihood function over this set; is the maximized 
log-likelihood. For Gaussian noise, the probability distri- 
bution function (pdf) of T is also particularly simple: a 
(perhaps non-central) distribution with 4 degrees of 
freedom. In their original paper by Jaranowski, Krolak 
& Schutz [l| (hereinafter referred to JKS), the J^-statistic 
was derived only for the case of a single GW detector and 
a smg le GW pulsar. Cutler & Schutz Q showed how the 
J^-statistic can be generalized in a straightforward man- 
ner to the cases of 1) a network of detectors noise curves, 
and 2) an entire collection of known sources. 

In practice, searches for nearly-monochromatic GWs 
arc separated into a few different types, depending on 
how much is know about the source. The different types 
of searches can have vastly different computational re- 
quirements. While the search for a GW counterpart 
to a known radio pulsar is trivial in terms of computa- 
tional burden, "all-sky" searches for GW pulsars with no 
known counterpart (and hence unknown frequency and 
frequency derivatives) are currently limited by the avail- 
able computational power. That is, we could dig deeper 



into the existing data sets if we possessed either larger 
computational resources or more efficient algorithms. In 
this paper we demonstrate a way of significantly improv- 
ing on the existing algorithms. 

Currently, the most sensitive all-sky searches are based 
on the following idea [J . For a G W pulsar with unknown 
frequency evolution (i.e., unknown /, /, etc.), computa- 
tional power required for an optimal (i.e, fully coherent, 
J^-bascd) search grows as a high-power of the total obser- 
vation time, T. Therefore in practice one divides T into 
some number N (typically of order 10^) short intervals 
of duration AT ~ T/N, performs an optimal, coherent 
search one each short interval, and then "adds up" the 
power from the all subintervals. More specifically, the 
current method is to calculate a "semi-coherent" detec- 
tion statistic J-'sc, defined as the sum of the T values from 
each of the short intervals: 

N 
i=l 

Because calcuting each Ti involves a maximization over 4 
free parameters, the pdf for 2Fsc is a distribution with 
4A^ d.o.f. But this is far more parameters than are actu- 
ally needed to describe the physical system! Consider the 
very first interval. The imbedded GW signal is described 
by 4 parameters: {Hq, But the triplet (/io,t, V') 

arc the same for all N intervals. All that changes from 
interval to interval is the overall phase 3>. Assuming 
maximum ignorance of the signal's phase evolution, one 
therefore needs iV — 1 additional phases to fully describe 
the signal. So the GW signal is fully described by only 
A'^ -t- 3 parameters. We define Tpr to be the maximized 
log- likelihood on this A^ -t- 3-dimensional space. Our main 
aims in this paper are (i) to demonstrate an efficient al- 
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gorithm for calculating Fpr and (ii) to and to illustrate 
its superiority as a detection statistic (superior in the 
sense of improved ROC curves, where ROC stands for 
"Receiver Operating Characteristics"). 

We note that our work is rather similar in spirtit to a 
recent paper by Dcrgachev Q , though we believe that we 
have advanced the idea considerably further. We note, 
too, a recent paper by Pletsch Q, which addresses the fol- 
lowng, rather different weakness in the detection statistic 
Eq. (jl.ip . To understand the problem, let (to7 ^i, • • • , t^) 
be the boundary points of the N short intervals, and con- 
sider any one such boundary time, say t^Q. In Eq. (|l.ip . 
data sampled ony slightly earlier than ^59 gets combined 
incoherently with data sampled only slightly later than 
^50, which clearly exacts a price in sensitivity. Pletsch 
overcomes this problem (which stems from the rather 
arbitrary choice of boundary points) using his "sliding 
window" technique We suspect that relatively sim- 
ple refinements of the detection statistic that we develop 
in this paper could also capitalize on Plctsch's basic in- 
sight, but we leave such refinements to future work. 

The plan of this paper is as follows. In Sec.|n]we briefly 
review the fully coherent J^-statistic, partly to establish 
notation. We generally try to align our notation with 
that of JKS, to ease comparison with their work. We also 
review the (currently used) semi-coherent J-"-statistic and 
its properties. 

As further motivation for our work, in Sec. II VI we con- 
sider a "warm-up" problem that we can easily treat ana- 
lytically and that qualitatively has much in common with 
our actual problem. In Sec.|Vl in order to illustrate both 
the use of J-pr, we present numerical results for one ex- 
ample search. For the same search, we also investigate 
the relative power/scnsitivty of J-pr versus J-ac- Our con- 
clusions are summarized in Sec. IVII 

This paper represents our "first-cut" analysis of the 
!Fpr statistic. There remains significant follow-up work to 
better elucidate the properties of J-pr and to implement 
it in realistic, hierarchical searches. This future work is 
also summarized in Sec. IVII 



II. REVIEW OF SIGNAL PROCESSING FOR 
GW PULSARS 

In this section we review the rudiments of signal pro- 
cessing that we will require, partly to fix notation. We 
also review both the coherent and semi-coherent versions 
of the J^-statistic. For simplicity of exposition, in this 
paper we will restrict the case of a single detector, and 
assume that the detector noise is stationary. The exten- 
sion to the more realistic case of multiple detectors with 
slowly changing noise spectra is completely straightfor- 
ward. 



A. Mathematics of signal processing 

We begin by reviewing the basic mathematics of sig- 
nal processing. For more details, we refer the reader to 
Thorne (1987), Finn & Chernoff (1993), and/or Cutler 
& Flanagan (1994) [IIIQ. 

Assuming that the noise is stationary and Gassuian, 
the noise spectral density Sh{f) determines a natural in- 
ner product (. . . I . . .) on the vector space of all detector 
outputs x{t): 

where x{f) and y(/) denote the Fourier transforms olx{t) 
and y(i), and Sh{f) is the single- sided spectral density of 
the noise. In terms of this inner product, the probability 
distribution function (pdf) for the noise n(i) takes the 
form 

pdf[n] = AAe-("l")/2, (2.2) 

where A/" is a normalization constant. Using ( ■ ■ ■ ) to 
denote " expectation value" (over many realizations of the 
noise), it follows from Eq. (|2.2p that 

((x|n)(y|n)) = (x|y). (2.3) 

In this paper, we will be concerned with waveforms 
h{t) that are nearly monochromatic (here meaning that 
their frequencies f{t) are slowly varying). In this case 
their inner product is equally simple in the time domain. 
Taking the measurement time interval to be to T, we 
have 

B. The fully coherent J^-statistic 

Next we briefly review the use of the coherent F- 
statistic in GW pulsar searches. For more details we 
refer the reader to Cutler & Schutz (2005) [1]. Consider 
a nearly monochromatic GW signal from an individual 
source with known sky location and known frequency 
evolution f{t). The GW signal is then characterized by 
four remaining unknowns: an overall amplitude A (equiv- 
alent to the combination /iq sin C,s\t?9 in the notation of 
JKS), two angles t and V' that characterize the waves' 
polarization (equivalent to determining the direction of 
the NS's spin axis), and an overall phase $. 

The GW signal h{t) registered by the detcctgor de- 
pends nonlinearly on t, $, but, crucially, one can make 
a simple change of variables-to (A^, A^, A^, A'') -such that 
dependence of hit) is linear in these new variables: 

4 

/i(t) = ^A"/ia(i) (2.5) 

a=l 
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where the four basis waveforms h"{t) arc defined by 

hi{t) = F+(t)cos$(t), h2(t) = Fx(t)cos$(t), 

h-sit) = F+(i)sin$(t), h4(t) = Fx (t)sin<I>(t) . (2.6) 

Here <f>(t) is the waveform phase at the deteetor: 



$(i) « 27r / fit')dt' 



(2.7) 



where f{t') is the measured GW frequency at the de- 
tector at time t' . The measured frequency includes the 
Doppler effect from the detector's motion relative to the 
source, as well as Einstein and Shapiro delays associated 
with the Earth's orbit around the Sun. When the GW 
pulsar is in a binary, then f(t') also includes the Roemer, 
Einstein, and Shapiro delays associated with that binary 
orbit. ( We emphasize that the known-pulsar searches 
described here do not require that the GW pulsar be iso- 
lated, but just that there exists an accurate timing model 
for the emitted waves.) The F+{t) and Fx{t) terms in 
Eq. (|2.6p are the beam-pattern functions describe the 
detector's response to the -I- and x polarizations, re- 
spectively. We note that the exact form of F+{t) and 
Fx (t) depends on one's convention for decomposing the 
waveform into "plus" and "cross" polarizations; a one- 
parameter family of choices is possible, corresponding to 
the freedom to rotate the axes around the line of sight. 
JKS follow the conventions of Bonazzola & Gourgoul- 
hon [13. 

Next we define the 4x4 matrix Tab by 



,dh_ dh_ 



= (ha I hb) 



(2., 



Because both the observation time and 1 day (the 
timescale on which the F" ^ (t) vary) are vastly 
larger than the period of the sought-for GWs (typ- 
ically 10^^ — lO^'^ s), we can replace eos^$(t), 
sin^$(t), and cos$(t)sin$(t) by their time-averages: 
cos^$(t), sin2$(t) i, while cos$(t)sin$(t) -J> 0. Then 
we have 



Til 

ri2 

^22 



F+{t) F+it)S^\fit))dt 
F+it) Fxit)S^\fit))dt 
Fx{t)Fxit)S^'{fit))dt; (2.9) 



additionally, T^^ « En, w ri2, « r22, and Fia « 

ri4 w F23 « r24 « 0. 

The best-fit values of A° satisfy 

^(x-^AXk-^A^hc) = (2.10) 



implying 



A° = E(r-'r''(x|hb) 



(2.11) 



Then 2T, which is defined to be twice the log of the 
maximized likelihood ratio, is just 

2J- = (x|x) - (x - ^ A'-hb I X - ^ A^e) 

6 c 

- E(r"'r'(x|ha)(x|hd). (2.12) 



Using 2J^ as one's detection statistic satisfies the 
Neyman-Pcarson criterion for an optimum test: it mini- 
mizes the false dismissal (ED) probability for any given 
false alarm (FA) probability. 

Writing x = n + h, and plugging into Eq. (|2.12p . we 
find 



(2^)=4+(h I h), 



(2.13) 



where we have used Eq. (|2.3p and the fact that 
((h I n)) = 0. More generally, it is easy to show that 
y = 2T follows a distribution with 4 degrees of free- 
dom (d.o.f) and non-centrality parameter = (h|h): 



^(2/) = x'(2/|4;p') 



(2.14) 



As pointed out by JKS, if we use the following com- 
plexified variables, 

2F, = (x|hi - zha) , 2Fb = (xjha - , (2.15) 



then the expression (|2.12[) for 2T can be re-written in a 
particularly simple form: Eq. (|2.12|) becomes 

2T ^^^\Faf + A\Fbf -2C^{FaFt)\ . (2.16) 
where 

A=(hi|hi), B = (h2|h2), C=(hi|h2), (2.17) 

and D = AB - C^. (Note that the A,B,C terms de- 
fined here are, in the single-detector case, larger than 
the A,B,C terms in JKS by a factor of the observation 
time T.) 

C. The "semi-coherent" 7^-statistic 

As mentioned above, the current method of incoher- 
ently combining the coherent results from successive in- 
tervals is just to sum of the J^-statistics from all the in- 
tervals: 



N 



2J~ sc — ^ ^ 2^ I 



(2.18) 



i=l 



It also easy to show that y = 2Tsc follows a distribu- 
tion with \N degrees of freedom: 



P(y)=x'(y|4iV;p2), 



(2.19) 
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where the non-centrahty parameter = Y^f^i pf- 

In the cases of interest to us, AN will generally be 
large, and then the distribution with AN d.o.f. can 
often be approximated as a Gaussian. Let y = 2J^, and 

let plt^HtiPi- Then 



P{v)=x\vm-pl,)^{^TTN) 



'i/2^-{y-<y>f/{SN) 



(2.20) 

where < y >= AN -\-p1^^. For example, using this approx- 
imation (and the fact that for a Gaussian P(?/), events 
2.326(7 above the mean occur 1% of the time), we see that 
the threshhold value yth that yields a 1% FA probability 



IS 



yth « 47V + 2.326> 



(large N) , 



(2.21) 



III. THE "PHASE-RELAXED" J^-STATISTIC 



We are now ready to define Fpr- Basically, Fpr co- 
incides with the full matched-filtering SNR^, under the 
assumption that the manifold of waveforms is + 3- 
dimensional (i.e., 4 parameters for the first segment, and 
TV — 1 for the relative phase offsets of the remaining seg- 
ments). What makes useful in practice is that we 
have also found a simple and efhcient method for calcu- 
lating it. 



A. Motivation and definition 

We begin by defining complex basis functions H-^- and 
by 



hi — i/i3 



(3.1) 



This complex representation is especially convenient 
for our purposes because and both transform very 
simply under an overall phase shift in $(i): under $(t) — > 
$(t) + 6, H+ and transform as H+{t) e-'^H+{t) 
and Hx{t) — )- e~"'ifx(i)- (Note that the minus sign in 
the exponent in the term e"*'' stems from the minus signs 
in the definitions of of iJ+ and Hy in Ea.f l3.ip . For these 
complexified signals, our usual inner product becomes a 
Hcrmitian one; for nearly monochromatic signals near 
frequency /, this Hcrmitian inner product is given simply 
by 



x|y = 



ShXf) 



x*{t)y{t)dt 



Clearly (x | y) = (y | x)*. 
Next we define Tap by 



= {Ha I Hp 



Pa 1 



(3.2) 



(3.3) 



where a and j3 run over -|-, x . It follows immediately that 
for GW data x{t), (twice the) J^-statistic is given by 



(r-i)"^(if„|x)(x|if,) 



(3.4) 



Now imagine breaking up the full integration time T 
into N intervals of duration AT^, for i = 1,2,--- ,A^. 
(We expect that in practice the AT^ will generally be of 
approximately the same length, but this is not required.) 
Next define x'^{t) to be the restriction of x{t) to the i*'' 
interval; i.e, x''{t) = x{t) for t in the i*'' interval, and 
x'^{t) = for t outside the i*'' interval. Then clearly we 
have 



2J" 



, aP 



(3.5) 



To motivate our definition of J-pr, recall that if we had 
practically limitless computer power at our disposal, then 
the most sensitive search would be a coherent matched- 
filter search over a fine grid covering the entire GW- 
pulsar parameter space. However for "blind" GW pul- 
sar searches (i.e., searches for GW pulsars whose sky 
location and/or time-changing frequency arc unknown), 
maintaining phase coherence between the template signal 
and true imbedded signal, over timescales of months to 
years, would require an extremely fine grid on parameter 
space, and (one easily shows) many of orders of magni- 
tude more computing power than is realistic [ll| 

The basic idea behind "semi-coherent" searches is to 
employ a detection statistic that is less sensitive to phase 
decohcrencc across the whole observation time, which 
allows one to use a much coarser grid on parameter 
space. In effect, one sacrifices some sensitivity in the 
interest of computational practicality. For our phase- 
relaxed J^-statistic, the idea is that the search-template 
signal should remain approximately in phase with the 
true, imbedded signal in each interval AT^-up to some 
constant phase "offset (J^-but that the 5i should be al- 
lowed to vary from interval to interval. That is, we re- 
place 

(r-^)"'(E.(^.l xO)(E, (xJ|i/,) 

{T-r'{i:. {Ha I x>^^-) (xJ |i/,)e-»*-^)3.6) 

Finally, we define (twice) Tpr to be the rhs of p.6p , max- 
imized over all phase-offsets Sii 



IT 



Si. 



(3.7) 

While there are N phase angles 5i , only — 1 of them 
are actually independent; i.e., it is easy to check that J-pr 
is invariant under Si ^ Si + c, where c is any constant. 



B. Maximizing over the phase offsets 5i 

The whole point of developing alternatives to the fully 
coherent J^-statistic is to save on computational cost, so 
for our phase-relaxed J^-statistic to be useful, we need 
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a reasonably efficient way of maximizing over tlie Si. In 
this section we demonstrate one efficient metliod. We 
demonstrate only the simplest version of this method, 
which we regard as basically an "existence proof" that 
efficient methods do exist. It should be clear by the end 
of this section that there are many variations on our ba- 
sic method by which one might attempt to improve its 
efficiency, but we defer such improvements to later work. 

Our method is as follows. We can simplify the appear- 
ance of the equations by defining 

K^'={r-Y\^'\Hp){H^\^') (3.8) 

and defining v to be the following A^-dimensional vector 
formed out of the phase offsets: 

VEE (e'■^^e"'^•••e"^"). (3.9) 

Note that K^' is Hermitian (i.e., K'^ = K^"), and that 
we can now re- write Eq. p.6|) as 

2Tr,r ^ max v.*K^'v^ . (3.10) 

Of course, K^^ is completely determined by the two com- 
plex templates and their inner products with the 
data, while our goal is to find the w,; that maximize 
VjK^'^Vi, subject to the N constraints that ViVi* — IVi. 
(We emphasize that i is not summed over in these con- 
straints.) Put another way: v lies on the unit N-torus 
(i.e, the unit circle cross itself N times.) Naturally, we 
employ the method of Lagrange multipliers to maximize 
v*K^^Vi on this constraint surface. Since there are N con- 
straints, we obtain N equations with N (real) Langrange 
multipliers Xj: 

K^\^ ^ XjVj Vj . (3.11) 

We emphasize that Eq. (|3.1ip is not an eigenvalue equa- 
tion, since in general the N values Xj will all be different. 

Next we find it convenient to introduce a pro- 
jection operator P operating on C^. Let w = 
(cie*^i , 026*''^ , • • • ,CAre'*"), where the Ci are all real. 
Then P is defined by 

Pw= (e'*Se*^%--- ,e*'^"). (3.12) 

I.e., the operator P takes any vector in and projects 
it down onto the unit torus. (Note that P is not a linear 
operator, but it is true that = P.) Then Eq. (|3.1ip is 
clearly equivalent to the requirement that 

PKw^w. (3.13) 

Hence the solution v is a fixed point of the operator 
PK. In fact, numerical experience shows that it is an 
attractive fixed point. That is, let Vq be some initial 
guess, and then operate on it repeatedly with PK. Define 
{PK^ = {PK){PK), {PKf EE {PK){PK){PK), etc. 
Then for Vq sufficiently close to the true solution vq, we 
find that 

{PKy vo V (3.14) 



as m increases. In practice, we find that the convergence 
is quite rapid, and that the initial guess Vq need not 
be particulary close to the solution vq. In numerical 
experiments (in many thousands of cases, and covering 
a large range of N) we found that the following initial 
guess always led to converge of the iterated sequence. 
For each segment ATi, it trivial to calculate the fully 
coherent J-i and the corresponding best-fit parameters 
for that segment alone: {Ai, Li,ipi^i). Then we take as 
our initial guess 

vo = (e'*\e**%---e**"); (3.15) 

i.e., the initial guess for the phase offset in each segment 
is the best-fit offset for that segment by itself. 

IV. ANALYTIC RESULTS FOR A RELATED, 
WARM-UP PROBLEM 

It is common sense that when one goes to solve some 
problem numerically, it is useful to have analytical results 
with which to compare it-ideally for a special case of the 
true problem, or, failing that, for some qualitatively sim- 
ilar problem. In this section we derive analytic results 
for the following case: Consider a vector space of wave- 
forms that is completely described by 2 parameters per 
interval-so 2N parameters in all, where N is large-and 
consider two different searches: one search that maxi- 
mizes the fit over those 2N parameters, and another, less 
efficient search, that begins with a 4A'^-dimensional vec- 
tor space (in which the true, 2iV-dim vector space lies), 
whose detection statistic is the maximized log-likelihood 
on the 4A^-dimensional space. That is, our two detection 
statistics are the 2N- and 4A^-dimensional J^-statistics, 
which in this section we will denote J-2n and J-^n- 

For each search, there is a threshold value pth such that 
the signal is detectable with FA = 0.01 and FD = 0.5. 
We can solve both problems at the same time, by consid- 
ering the general M-dimensional search. Then the expec- 
tation value of J^M is {^m) = M + and its standard 
deviation is ctm = {2M + 4p^)^/^. For large M, the 
function approaches a Gaussian, so we will approximate 
the pdf of J^M as a Gaussian with this mean and stan- 
dard deviation. Then the threshold for detection with 
FA = 0.01 is 

Pm ^ {J'M)+V2aM CTk~\2FA) (4.1) 

= M + 3.29M^^^ , (4.2) 

where in Eq. (??) both J-i\j and ctm are to be evaluated 
at p = 0. Therefore pf^ = V3.29AF/^ = l.SUM^/'^, and 
we have 

pI%/pI% = 2'/'^1.189. (4.3) 

Therefore using the correct statistic allows one to see 
sources 19% farther away. 

For comparison with results in the next section, we also 
plot in Fig. 1 the FA vs. FD curves for the two statistics, 
for a range of p values. 
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log^ijFA vs. FD for different values ofp 




FD 

FIG. 1: Compares the false alarm (FA) probabilities for the 
two detection statistics, T4N and J-2n, as a function of false 
dismissal (FD) probability, for several values of the squared 
signal strength, p^. The blue curves are for J-4jv and the green 
for J-2N- From upper to lower, the squared signal strengths 
are = 25, 50, 75, 100. 



V. NUMERICAL RESULTS FOR ONE 
EXAMPLE SEARCH 

To illustrate the utility of our Jp,. statistic, in this 
section we present results for a simple, one-parameter 
family of examples (where the varied parameter is the 
strength of the embedded GW signal) , and we compare 
the the effectiveness of J-pr and J-gc- 

We fix the number of intervals at = 100, and eval- 
uate search effectiveness for signals with a range of total 

= Si^i Pf • ■^ill eventually consider a range of p, 
but for now imagine p as fixed. For simplicity, in this ex- 
ample we will consider a case where the pi are the same 
for all i, so pf = p^ /N , and where the T^j^ matrices are 
also the same for all i: = 3, T^^^ = 1 = rx_|_ and 

r*xx = 1 for all i. 

We decompose the measured signal into waveform 
plus noise, 

x'^h' + n\ (5.1) 

For each filtering the data with iJ+ and Hx pro- 
duces two complex numbers: c!|_ = (x' | and c'^ = 
(x' I iJx ) • Clearly, the measured signals can be de- 
composed as 

= (h' I + (n- 1 (5.2) 

^ .9L + <- (5.3) 

We simulate the noise piece by taking random draws 
of (pairs of) complex numbers from a Gaussian distribu- 



tion with covariance matrix 

^ ((i/Jn-) <{n^\Hp)) 

= (i/„ \H0)^T^p (5.4) 

Again for simplicity, we will consider a case where the gf 
are the same for each i, modulo a random, complex phase 
factor. Our particular (and rather arbitrary) choice is 

[ .9+ , gf] = {2N)-'l^p [ 2 + V6 , V6] e^^" , (5.5) 

where the are random phases drawn uniformly from 
[0,27r). One easily checks that E^iC*' I ^0 = ■ The 
inclusion of the e*"^' terms reflects our goal of model- 
ing a case where frequency evolutions of the tempate 
and the true signal are so mismatched that their rela- 
tive phases jump significantly and randomly from one 
interval to the next. Choosing the (p* randomly cor- 
responds to the " worst-case scenario" , where the true- 
versus-template phase offsets show no pattern. In prac- 
tice, we expect that the situation will often be much more 
favorable for searches: i.e., the phase offsets might very 
often be well fit by some low-order polynomial in time. In 
a later paper we plan to investigate the extent to which 
the time-evolution of the offsets can be fit by a few pa- 
rameters, and how that information can be exploited to 
speed up other parts of the search. 

Given one simulated data set (200 complex num- 
bers), we compute 2Tpr. We repeat for 10000 data sets 
to determine the distribution of 2Tpr, and calculate its 
mean (2J-pr{p))^ standard deviation apr{p), skewness, 
and kurtosis. In practice, we find that the skewness and 
kurtosis are relatively small (as might be expected, since 
our N is large), so for the rest of this section we will 
approximate p{2J^pr; p^) as simply a Gaussian with our 
measured mean and standard deviation. Given these dis- 
tributions it is completely straightforward to determine 
the false alarm probability FA for any threshold value 
2Tp!^, and to calculate the false dismissal probability FD 
for any pair of 2 and p^ . 

We expect that, in practical searches, J^^Jf will find its 
main use in hierarchical search algorithms, in which a 
very coarse search at relatively low threshold identifies 
candidates for further examination, and these are win- 
nowed down in successive stages [H, [l^. In this context, 
one generally wants a fairly small FD rate (< 1%, say), 
so as not to lose any events, and strongly prefers a very 
low FA probability, to reduce the computational cost of 
follow-ups. With this application in mind, in Figs. 2 and 
3. wc plot FA as a function of FD, for several values of 

How much does employing Tpr increase the sensitity 
of a blind GW pulsar search? Obtaining a useful and 
accurate answer to this question is much more compli- 
cated than it might initially seem, since the most sen- 
sitive known search algorithms for GW pulsars are hi- 
erarchical searches [l^, [13 ■ These searches involve in- 
volve several stages, with successive stages "ruling out" 
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log^QpA vs. FD for different values of p 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 
FD 



FIG. 2: Compares the false alarm (FA) probabilites for the 
two detection statistics, J^sc and Tpr, as a function of false 
dismissal (FD) probability, for several values of the (square 
of) the signal strength, . The blue curves are for Tsc and 
the green for J-pr- From upper to lower, the signal strengths 
are = 75, 100, 125. 



log^ijFA vs. FD for different values of P 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 
FD 



FIG. 3: Same as in Fig. 2, except that, from the upper to 
lower curves, the signal strengths are = 150, 175, 200. 



an ever-increasing fraction of the parameter space. We 
imagine that the most sensitive search-at some fixed, re- 
alistic computational cost-might use J-pr for only some 
of its stages. And while calculating J-pr clearly requires 
more floating point operations than calculating J"sc, it is 
premature to compare these costs in a detailed way, since 
(i) there has as yet been no attempt to speed up our it- 
erative relaxation scheme, and (ii) the extra information 
that comes with the phase offsets 5i can presumably be 



used to speed up other parts of the search. Despite these 
difficulties we can obtain a rough estimate of the sensi- 
tivity improvement from using p^^, as follows. We have 
seen that given some detection statistic S, one naturally 
obtains a map from the total to curves in the FA—FD 
plane. Let C{S] p^) denote that curve. Then we can look 
for pairs p^^ and p1^ such that C{Fpr\ p^r) ^^^^ close to 
C{Fsc] pIc)- Three such pairs are shown in Fig. 4. We 
see that, in the most relevant portion of FD — FA plane, 
C{Fpr] 140) lies close to C{Fsc] 170), C{Fpr; 160) lies close 
to C(J"sc; 200), and C{Tpr] 190) lies close to C{Fsc\ 250). 
Thus, based on this example, we might estimate that 
using J^pr rather than affords an increase in sensitiv- 
ity of - 20 - 30% in or - 10 - 15% in p. Clearly, 
to obtain a more reliable estimate we should perform a 
Monte Carlo simulation (based on random locations of 
the source on the sky, and random orientations of the 
GW pulsar's spin axis). We plan to do this in follow-up 
work. 

Comparing Fig. 1 to Fig. 2, we see that the sensitiv- 
ity gain from replacing the detection statistic Tsc with 
Fpr is qualitatively similar to the gain from J\4jv J'2N, 
but that the latter gain is greater (at least for our one 
numerical example, and for total SNR ~ 10). That may 
seem surprising, since in the former case we are elimi- 
nating 3iV — 3 redundant parameters, while in the latter 
case we are eliminating only 2A^ redundant parameters. 
We conjecture that the main reason that the replacement 
J-'sc — >■ J^pr "buys us less" in sensitivity is the following. 
In calculating J-2n, the noise contributions from differ- 
ent intervals i still get combined incoherently. However 
in J-pri the maximization over the phase offsets di allows 
the noise contributions to combine coherently. Indeed for 
Pi <^ 1, the offsets 6i are determined more by the noise 
than by the imbedded waveform. However we leave a 
thorough investigation of this effect to future work. 



VI. SUMMARY AND FUTURE WORK 



In this paper we defined a "phase-relaxed" J^-statistic, 
denoted J-pr , to be used in cases where the total observa- 
tion time is sufficiently long that straightforward calcu- 
lation of the fully coherent J^-statistic over the relevant 
parameter space is computationally intractable. The cal- 
culation of J^pr takes as input the results from coherent 
searches over N shorter time intervals. Our Tpr coincides 
with the fully coherent J^-statistic under the approxima- 
tion that the phase offsets between template and imbed- 
ded signal are treated as an additional A'^ — 1 independent 
parameters. We also demonstrated one efficient, itera- 
tive method for calculating Tpr- We regard our iterative 
method as an "existence proof for efficient algorithms. 
In future work we intend to explore variations on our 
basic method that we suspect would lead to substantial 
improvements in computational cost. We illustrated the 
use of J-pr in one simple family of examples, in which the 
sensitivity improvement (compared to J^sc) was shown to 
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log^QpA vs. FD for different values of p 
-2 I 1 , , , 1 1 , — 




-14' ' • • • ' ' • • • 1 

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 
FD 

FIG. 4: Similar to Figs. 2 and 3, except that here the sensi- 
tivity of the Tsc statistic is compared to that of Tpr at lower 
p^. From upper to lower, the signal strengths for Tbc are for 
= 140, 160, 190. 



be ~ 10 - 15%. 

Our example was based on the "worst-case" assump- 
tion that where the phase offsets Si are completely ran- 
dom. In follow-on work we intend to examine the more 
realistic case where the 6i are a (reasonably low-order) 
polynomial in time, and we plan to calculate the in- 
creased sensitivity based on a very large number of cases. 



in Monte Carlo fashion. 

Other follow-on projects that we intend to work on 
include i) the development of new vetoes [l^ for instru- 
mental artifacts, since, e.g., for true GW pulsars the pa- 
rameters {ho, L, ip) calculated from the first half of the ob- 
servation should be consistent with those calculated from 
the second half; ii) the use of the Si to quickly converge 
on improved estimates for the GW pulsar's frequency and 
spindown parameters, and iii) the optimal use m 
multi-stage, hierarchical searches. 

Finally, while our primary interest in this paper has 
been the application of the J^pr-statistic to GW pulsars, 
it has not escaped our notice that the same idea and for- 
malism can be applied, with only trivial modifications, to 
searches for (quasi-circular, non-precessing) inspiraling 
binaries. With binaries, we expect J^pr to be most useful 
in those cases where the observed GW signal has a very 
large number of cycles; e.g., searches for neutron-star bi- 
naries by proposed GW detectors that have reasonable 
sensitivity at ^ 1 Hz, such as the Einstein Telescope, 
Decigo, or the Big Bang Observer. 
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